{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 多量子比特含噪模拟器\n",
    "*版权所有 (c) 2021 百度量子计算研究所，保留所有权利。*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 内容概要\n",
    "\n",
    "在本教程中，我们将介绍如何使用脉冲层面的多量子比特模拟器。本教程的概要如下：\n",
    "\n",
    "- 背景介绍\n",
    "- 准备工作\n",
    "- 含噪模拟器的门级控制\n",
    "- 含噪模拟器的脉冲级控制\n",
    "    - 构建系统的物理模型\n",
    "    - 拉比振荡\n",
    "    - Cross-Resonance 效应\n",
    "    - 通过 Ramsey 实验刻画 ZZ 串扰\n",
    "- 总结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 背景介绍\n",
    "\n",
    "在脉冲层面模拟量子比特的演化能让我们更深入地了解量子门的工作原理与噪声的来源。在超导电路中，transmon 量子比特通过施加微波或磁通脉冲来控制。然而，量子门的性能通常受到各种因素的影响——量子比特与环境相互作用导致的退相干、串扰噪声以及向高能级的泄漏。\n",
    "\n",
    "量脉提供的多比特含噪模拟器使我们能够在多量子比特的含噪声设备上模拟量子操作，以更好地了解量子计算背后的物理原理。该含噪模拟器包含了几种主要的噪声：退相干噪声、脉冲幅值失真噪声以及串扰噪声。基于该含噪模拟器我们将重点呈现超导量子计算测控流程中几项常见的操作，包括：拉比振荡实验，Cross-Resonance 效应以及通过 Ramsey 实验刻画 ZZ 串扰。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 准备工作\n",
    "\n",
    "成功安装量脉后，您可以按照本教程运行下面的量脉程序。要运行此教程，您需要从 Quanlse 和其它常用的 Python 库导入以下包："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "from Quanlse.remoteOptimizer import remoteOptimize1Qubit as optimize1q\n",
    "from Quanlse.remoteSimulator import remoteSimulatorRunHamiltonian as runHamiltonian\n",
    "from Quanlse.Superconduct.Simulator import PulseModel\n",
    "from Quanlse.Superconduct.Simulator.PulseSim3Q import pulseSim3Q\n",
    "\n",
    "from Quanlse.QWaveform import QJob, QJobList\n",
    "from Quanlse.QOperator import driveX, driveY, a, sigmaZ, number, driveZ\n",
    "from Quanlse.QWaveform import square, gaussian\n",
    "from Quanlse.Utils.Functions import basis, tensor, expect, dagger, partialTrace, project, computationalBasisList, population\n",
    "from Quanlse.Utils.Bloch import rho2Coordinate, plotBloch\n",
    "from Quanlse.Utils.Plot import plotBarGraph\n",
    "from Quanlse.QOperation.FixedGate import H, X, Y, CR, Z, CNOT\n",
    "from Quanlse.QOperation.RotationGate import RZ\n",
    "from Quanlse.Superconduct.SchedulerSupport import SchedulerSuperconduct\n",
    "from Quanlse.Utils.Infidelity import unitaryInfidelity\n",
    "from Quanlse.Superconduct.SchedulerSupport.PipelineCenterAligned import centerAligned\n",
    "\n",
    "from math import pi\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用量脉云服务之前，我们需要一个 token 来访问云端。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Import Define class and set the token for cloud service\n",
    "# Please visit http://quantum-hub.baidu.com\n",
    "from Quanlse import Define\n",
    "Define.hubToken = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 含噪模拟器的门级控制\n",
    "\n",
    "\n",
    "我们可以在量子门层面进行模拟。在这一部分，我们需要使用预先定义、具有默认配置的 `PulseModel()` 对象。\n",
    "\n",
    "为了创建一个三量子比特物理系统，我们首先通过调用 `pulseSim3Q()` 来实例化 `PulseModel()` 对象。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "model = pulseSim3Q(frameMode='lab', dt=0.01)\n",
    "model.savePulse = False\n",
    "model.pipeline.addPipelineJob(centerAligned)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了使用量脉调度器来定义一个产生 GHZ 态的量子电路， 我们可以通过调用 `gate(model.Q[index])` 将这些量子门添加到物理模型中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Hadamard gate \n",
    "H(model.Q[0])\n",
    "\n",
    "# CNOT\n",
    "CNOT(model.Q[0], model.Q[1])\n",
    "CNOT(model.Q[1], model.Q[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面定义的量子电路的脉冲序列是通过调用函数 `model.schedule` 生成的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "scheJob = model.schedule()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定义该三个量子比特初始态都在基态，也就是系统初始态为 $|\\psi\\rangle = |000\\rangle$，然后运行模拟，并且画出每个测量结果的概率分布。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Run the simulation\n",
    "res = model.simulate(job=scheJob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the result\n",
    "popList = [abs(item ** 2) for item in res[0]['state'].T[0]]\n",
    "basisList = computationalBasisList(3, 3)\n",
    "plotBarGraph(basisList, popList, \"Result\", \"Outcome\", \"Population\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上图可见，由于噪声的存在，测量结果会得到 $000$ 和 $111$ 外的值。使用含噪模拟器研究退相干情况下的模拟，用户可以尝试设置 `runHamiltonian()` 里的参数 `isOpen` 为 `True`，得到演化过程的密度矩阵，但是会花费一定的时间得到结果。如果想知道退相干如何影响超导量子计算，可以参考教程：[单量子比特含噪模拟器](https://quanlse.baidu.com/#/doc/tutorial-single-qubit-noisy-simulator)。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 含噪模拟器的脉冲级控制\n",
    "\n",
    "含噪模拟器支持脉冲级的模拟。通过对量子比特 $q_i$ 输入脉冲波形以及其它参数，就能得到系统的动力学演化，这使得我们可以从更底层模拟超导量子计算的控制环节。这里我们利用该模拟器，演示几个真实实验中常用的操作。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 构建系统的物理模型\n",
    "\n",
    "通常，我们可以用 Duffing 振子模型来描述超导量子比特的物理模型。以三比特模型为例，当量子比特 $q_0$ 与 $q_1$、$q_1$ 与 $q_2$ 之间有耦合的时候，在实验坐标系，哈密顿量可以表示为:\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\sum_{i=0}^2 \\omega_i \\hat{a}^\\dagger_i \\hat{a}_i + \\sum_{i=0}^2 \\frac{\\alpha_i}{2} \\hat{a}^\\dagger_i \\hat{a}^\\dagger_i \\hat{a}_i \\hat{a}_i + g_{01} (\\hat{a}^\\dagger_0 \\hat{a}_1 + \\hat{a}_0 \\hat{a}^\\dagger_1) +  g_{12} (\\hat{a}^\\dagger_1 \\hat{a}_2 + \\hat{a}_1 \\hat{a}^\\dagger_2),\n",
    "$$\n",
    "\n",
    "其中，参数 $\\omega_i$、$\\alpha_i$ 分别是量子比特 $q_i$ 的本征频率以及失谐量；而 $g_{ij}$ 则是量子比特 $q_i$ 以及量子比特 $q_j$ 之间的耦合强度。算符 $\\hat{a}_i$、$\\hat{a}^\\dagger_i$ 是量子比特 $q_i$ 的湮灭算符和产生算符。\n",
    "\n",
    "在本教程中，我们构造了一个包含三个量子比特的含噪模拟器。首先，我们定义一些硬件参数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "qubitNum = 3  # The number of qubits\n",
    "level = 3  # The energy level for each qubit\n",
    "\n",
    "anharm = -0.33 * 2 * pi  # The anharmonicity of the qubit, in 2 pi GHz\n",
    "wq0 = 4.914 * 2 * pi  # The frequency for qubit 0, in 2 pi GHz \n",
    "wq1 = 5.100 * 2 * pi  # The frequency for qubit 1, in 2 pi GHz\n",
    "wq2 = 5.200 * 2 * pi  # The frequency for qubit 2, in 2 pi GHz\n",
    "g01 = 0.0038 * 2 * pi  # The coupling strength of the interaction between qubit 0 and qubit 1, in 2 pi GHz\n",
    "g12 = 0.0020 * 2 * pi  # The coupling strength of the interaction between qubit 1 and qubit 2, in 2 pi GHz\n",
    "\n",
    "dt = 1.  # The sampling time of AWG\n",
    "\n",
    "# T1 relaxation time for qubit 0, qubit 1, and qubit 2, in nanoseconds\n",
    "t01 = 1000  \n",
    "t11 = 1120\n",
    "t21 = 1300\n",
    "\n",
    "# T2 dephasing time for qubit 0, qubit 1, and qubit 2, in nanoseconds\n",
    "t02 = 500\n",
    "t12 = 450\n",
    "t22 = 600\n",
    "\n",
    "# The random amplitude distortion\n",
    "ampNoise = 0.02"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如下，我们可以通过实例化 `PulseModel` 类中的一个对象来创建该物理系统。物理系统中的噪声包含 $T_1$ 弛豫噪声, $T_2$ 失相噪声与振幅噪声。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "qubitFreq = {0: wq0, 1: wq1, 2: wq2}  # Qubit frequency for each qubit\n",
    "qubitAnharm = {0: anharm, 1: anharm, 2: anharm}  # Qubit anharmonicity for each qubit\n",
    "qubitT1 = {0: t01, 1: t11, 2: t21}  # Relaxation time \n",
    "qubitT2 = {0: t02, 1: t12, 2: t22}  # Dephasing time\n",
    "couplingMap = {(0, 1): g01, (1, 2): g12}  # Coupling map\n",
    "\n",
    "# Create an instant of PulseModel\n",
    "model = PulseModel(subSysNum=qubitNum,\n",
    "                   sysLevel=level,\n",
    "                   qubitFreq=qubitFreq,\n",
    "                   qubitAnharm=qubitAnharm,\n",
    "                   couplingMap=couplingMap,\n",
    "                   T1=qubitT1,\n",
    "                   T2=qubitT2,\n",
    "                   dt=dt,\n",
    "                   ampSigma=ampNoise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "于是一个包含三个超导量子比特和三种常见噪声的模拟器就构建完成了。下一步，我们通过 `createQHamiltonian()` 方法创建一个 `QHamiltonian` 的对象。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "ham = model.createQHamiltonian()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross-Resonance 效应\n",
    "\n",
    "全微波控制是实现超导电路量子控制的方案之一。 在该方案中，双量子比特门的实现利用了两个弱耦合量子比特的 Cross-Resonance 效应。该效应通过在控制量子比特上施加以目标量子比特频率为驱动频率的脉冲完成。理想情况下，控制量子比特和目标量子比特之间的 $\\hat{\\sigma}_z \\otimes \\hat{\\sigma}_x$ 相互作用在所有类型的相互作用中占主导地位 \\[1\\]。如果想知道更多关于 CR 门的内容，可以参考教程：[Cross-Resonance 门](https://quanlse.baidu.com/#/doc/tutorial-cr)。\n",
    "\n",
    "在我们的模拟中，我们再次以不同的脉冲幅值驱动量子比特 $q_0$ （控制量子比特），脉冲的驱动频率为量子比特 1 的频率。这一过程可以由方法 `addWaveRot(index, waves, detuning)` 实现，其中 `index` 是作用的量子比特编号、`waves` 是脉冲波形、 `detuning` $\\Delta$ 是频率差（$\\Delta = \\omega_i - \\omega_d$，其中 $\\omega_i$ 是量子比特 $q_i$ 的频率而 $\\omega_d$ 是驱动频率）。 \n",
    "\n",
    "接下来我们改变脉冲幅值，观察在 cross-resonance 效应下不同幅值下每个量子比特的布居数变化。本例以量子比特 $q_1$ 为控制比特，$q_0$ 为目标比特，也就是说用驱动频率是 $q_0$ 的本征频率的脉冲驱动控制比特 $q_1$。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "dCoef = 0.03 * (2 * pi)  # The drive strength of the pulse\n",
    "ampCR = np.linspace(0, 0.5, 40)  # The amplitudes in arbitrary unit \n",
    "amps = ampCR * dCoef  \n",
    "detuning = wq1 - wq0  # The detuning of the pulse\n",
    "\n",
    "# jobList = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt, title='cr')\n",
    "jobList = ham.createJobList()\n",
    "\n",
    "# Fix the gate time\n",
    "tg = 950\n",
    "\n",
    "# Append each job to the jobList\n",
    "for amp in amps:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job = ham.createJob()\n",
    "    job.addWaveRot(1, waves=square(tg, amp), t0=0., detuning=detuning)  # Apply pulse at qubit 1\n",
    "    job = model.getSimJob(job)\n",
    "    jobList.addJob(jobs=job)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在，我们设置该三比特系统的初始态为 $|\\psi\\rangle = |010\\rangle$，即控制比特 $q_1$ 处在激发态，然后进行模拟。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Define the initial state of |010>\n",
    "psi0 = tensor(basis(level, 0), basis(level, 1), basis(level, 0))  \n",
    "\n",
    "# Run the simulation\n",
    "result = runHamiltonian(ham=ham, state0=psi0, jobList=jobList)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定义每个量子比特 $q_i$ 上的第一激发态的投影算符。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "prj01 = tensor(basis(3, 1) @ dagger(basis(3,1)), np.identity(level), np.identity(level))  \n",
    "prj11 = tensor(np.identity(level), basis(3, 1) @ dagger(basis(3,1)), np.identity(level))  \n",
    "prj21 = tensor(np.identity(level), np.identity(level), basis(3, 1) @ dagger(basis(3,1)))  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算每个量子比特上投影算符 $|1\\rangle\\langle1|$ 的期望值，然后绘制激发态布居数随脉冲幅值变化的图。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Initialize the list of expected values\n",
    "num0List = []\n",
    "num1List = []\n",
    "num2List = []\n",
    "\n",
    "for res in result.result:\n",
    "    state = res['state']  # The final state of each job\n",
    "    num0Expect = expect(prj01, state)  # Compute the expected values of the projector |1><1|\n",
    "    num1Expect = expect(prj11, state)\n",
    "    num2Expect = expect(prj21, state)\n",
    "    num0List.append(num0Expect)\n",
    "    num1List.append(num1Expect)\n",
    "    num2List.append(num2Expect)\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "# Plot the figure of CR effect\n",
    "plt.plot(ampCR, num0List, label='qubit0')\n",
    "plt.plot(ampCR, num1List, label='qubit1')\n",
    "plt.plot(ampCR, num2List, label='qubit2')\n",
    "\n",
    "plt.xlabel('Amplitudes (a.u.)')\n",
    "plt.ylabel(r'Population of $|1\\rangle$')\n",
    "plt.title('Cross-Resonance effect')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可见，目标比特 $q_0$ 的激发态布居数随着脉冲幅度增加而改变，而控制比特 $q_1$ 在脉冲幅度较小的时候基本处于激发态，与 $q_1$ 耦合的比特 $q_2$ 则一直处于基态。当脉冲幅度进一步增大的时候，也会不可避免地影响到到量子比特 $q_1$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 通过 Ramsey 实验刻画 ZZ 串扰\n",
    "\n",
    "ZZ 串扰是耦合量子比特间寄生相互作用的主要来源，它源于量子态更高能级间的相互作用。对于一对直接或者间接耦合的量子比特，其计算空间的有效哈密顿量为 \\[2\\]：\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\rm eff} = \\omega_{0}\\frac{\\hat{\\sigma}_{z}^0 \\otimes I_1}{2} + \\omega_{1}\\frac{I_0\\otimes\\hat{\\sigma}_{z}^1}{2} + \\xi \\frac{\\hat{\\sigma}_{z}^0 \\otimes \\hat{\\sigma}_{z}^1}{2},\n",
    "$$\n",
    "\n",
    "其中 $\\omega_0$，$\\omega_1$ 是量子比特的频率，$\\xi$ 是 ZZ 串扰的强度。$\\xi$ 由量子态 $|11\\rangle \\leftrightarrow |10\\rangle$ 与 $|01\\rangle \\leftrightarrow |00\\rangle$ 间的跃迁频率定义:\n",
    "\n",
    "$$\n",
    "\\xi = \\left(E_{11} - E_{10}\\right) - \\left(E_{01} - E_{00}\\right),\n",
    "$$\n",
    "\n",
    "其中，$E_{ij}$ 是态 $|ij\\rangle$ 的能级。实际上，我们可以通过 Ramsey 实验探测到这个频率差，因此对该串扰进行度量与刻画。在逻辑量子门层面，这相当于在两个 $H$ 门间有一段空转的时间 \\[3\\]。\n",
    "\n",
    "为了更好地描述 ZZ 串扰带来的影响，我们定义一个耦合强度更大的三量子比特模型（耦合强度约为 6 ~ 40 MHz）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "dt = 0.2  # The sampling time\n",
    "level = 3  # The system level\n",
    "qubitNum = 3  # The number of qubits\n",
    "\n",
    "g01 = 0.0377 * (2 * pi)\n",
    "g12 = 0.0060 * (2 * pi)\n",
    "\n",
    "# Coupling map\n",
    "couplingMap = {\n",
    "    (0, 1): g01,\n",
    "    (1, 2): g12\n",
    "}\n",
    "\n",
    "# Qubits frequency anharmonicity\n",
    "anharm = - 0.33 * (2 * pi)\n",
    "qubitAnharm = {0: anharm, 1: anharm, 2: anharm}  # The anharmonicities for each qubit\n",
    "\n",
    "# Qubit Frequency\n",
    "qubitFreq = {\n",
    "            0: 5.5904 * (2 * pi),\n",
    "            1: 4.7354 * (2 * pi),\n",
    "            2: 4.8524 * (2 * pi)\n",
    "            }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过类 `PulseModel()` 创建物理模型，并且根据模型创建哈密顿量 `ham`。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "model = PulseModel(subSysNum=qubitNum, sysLevel=level, couplingMap=couplingMap,\n",
    "                    qubitFreq=qubitFreq, dt=dt, qubitAnharm=qubitAnharm)\n",
    "    \n",
    "ham = model.createQHamiltonian()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过调用单量子比特门的优化器, 生成作用在不同量子比特的 $H$ 门与 $X$ 门。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Define function to generate the QJob for gate of specified qubit \n",
    "def generateGate(gate, index):\n",
    "    job1q, _ = optimize1q(ham=ham.subSystem(index), uGoal=gate.getMatrix(), targetInfid=1e-5)\n",
    "    job3q = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    waves = job1q.waves\n",
    "    ops = job1q.ctrlOperators\n",
    "    for key, op in ops.items():\n",
    "        job3q.addWave(operators=op, onSubSys=index, waves=waves[key])\n",
    "     \n",
    "    return job3q\n",
    "\n",
    "# Generate the gates needed\n",
    "h0 = generateGate(H, 0)  # H gate on qubit 0 \n",
    "h1 = generateGate(H, 1)  # H gate on qubit 1\n",
    "x1 = generateGate(X, 1)  # X gate on qubit 1\n",
    "x2 = generateGate(X, 2)  # X gate on qubit 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "maxTime = 500  # The delayed time in Ramsey experiment, in nanosecond.\n",
    "freq = 3 / maxTime  # Detuning. \n",
    "\n",
    "# Generate job for delayed time \n",
    "def generateIdle(tg, index):\n",
    "    jobIdle = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    jobIdle.appendWave(operators=driveZ, onSubSys=index, waves=square(tg, 2 * pi * freq))\n",
    "    return jobIdle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定义两个不同的 `jobList` 对象——一个初始量子态是 $|00\\rangle$，另一个通过作用一个 $X$ 门到量子比特 $q_1$ 上将其初始态是 $|01\\rangle$。然后在量子比特 $q_0$ 上进行 Ramsey 实验。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# jobList with initial state |00>\n",
    "jobListGrd = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "\n",
    "# jobList with initial state |01> (by applying X gate) \n",
    "jobListExd = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "\n",
    "# Define the delayed time\n",
    "tgList = np.linspace(0, maxTime, 50)\n",
    "\n",
    "# Define jobList with initial state |00>\n",
    "for tg in tgList:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job += h0\n",
    "    job += generateIdle(tg, 0)\n",
    "    job += h0\n",
    "    jobListGrd.addJob(job)\n",
    "\n",
    "# Define jobList with initial state |01>\n",
    "for tg in tgList:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job += x1\n",
    "    job += h0\n",
    "    job += generateIdle(tg, 0)\n",
    "    job += h0\n",
    "    jobListExd.addJob(job)\n",
    "\n",
    "# Run the simulation\n",
    "stateInit = tensor(basis(level, 0), basis(level, 0), basis(level, 0))\n",
    "resultGrd = runHamiltonian(ham, state0=stateInit, jobList=jobListGrd)\n",
    "resultExd = runHamiltonian(ham, state0=stateInit, jobList=jobListExd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "绘制量子比特 $q_0$ 的第一激发态的布居数随空转时间变化的分布。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "num0List = []\n",
    "num1List = []\n",
    "\n",
    "# projector |1><1| of qubit 0\n",
    "prj1 = tensor(basis(level, 1) @ dagger(basis(level, 1)), np.identity(9))\n",
    "\n",
    "# append the result to the list\n",
    "for res0, res1 in zip(resultGrd, resultExd):\n",
    "    psi0, psi1 = res0['state'], res1['state']\n",
    "    rho0, rho1 = psi0 @ dagger(psi0), psi1 @ dagger(psi1)\n",
    "    num0List.append(expect(prj1, rho0))\n",
    "    num1List.append(expect(prj1, rho1))\n",
    "\n",
    "\n",
    "# plot the result\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(tgList, num0List, '.b', label=r'$|00\\rangle$')\n",
    "plt.plot(tgList, num1List, '.r', label=r'$|01\\rangle$')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel('Population of excited state of qubit 0')\n",
    "plt.title('Ramsey on Q0')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ZZ 串扰的强度可以通过不同的 Ramsey 振荡频率测量。因此，我们通过用余弦函数去拟合得到的结果，并且计算得到频率 $f_1$，$f_2$。串扰强度为 $\\xi / \\left( 2\\pi \\right) = |f_1 - f_2|$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Define the fitting curve\n",
    "def fit(x, omega, theta):\n",
    "    return - 0.5 * np.cos(omega * x + theta) + 0.5\n",
    "\n",
    "# Fit the curve\n",
    "para1Fit, _ = curve_fit(fit, tgList, num0List, [2.1 * pi * freq, 0])\n",
    "para2Fit, _ = curve_fit(fit, tgList, num1List, [2 * pi * freq, 0])\n",
    "step = 0.01\n",
    "y1Fit = [fit(x, para1Fit[0], para1Fit[1]) for x in np.arange(tgList[0], tgList[-1], step)]\n",
    "y2Fit = [fit(x, para2Fit[0], para2Fit[1]) for x in np.arange(tgList[0], tgList[-1], step)]\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "# Plot the curve\n",
    "plt.plot(np.arange(tgList[0], tgList[-1], step), y1Fit)\n",
    "plt.plot(np.arange(tgList[0], tgList[-1], step), y2Fit)\n",
    "plt.plot(tgList, num0List, '.b', label=r'$|00\\rangle$')\n",
    "plt.plot(tgList, num1List, '.r', label=r'$|01\\rangle$')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel('Population of excited state of qubit 0')\n",
    "plt.title('Ramsey on Q0')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the crosstalk strength\n",
    "xiEst = abs(para1Fit[0] - para2Fit[0]) \n",
    "print(f'Coupling strength: {g01 * 1e3 / (2 * pi)} MHz')\n",
    "print(f'ZZ crosstalk strength: {xiEst * 1e3 / (2 * pi)} MHz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于量子比特 $q_0$ 以及 $q_1$ 之间存在比较强的耦合强度，能观察到频率差比较大，也就是 ZZ 串扰比较大。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以验证较小的耦合强度下的 ZZ 串扰也会较小。重复模拟该实验以得到量子比特 $q_1$ 与 $q_2$ 间的 ZZ 串扰强度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# jobList with initial state |00>\n",
    "jobListGrd = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "\n",
    "# jobList with initial state |01> (by applying X gate)\n",
    "jobListExd = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "\n",
    "# Define the delayed time\n",
    "tgList = np.linspace(0, maxTime, 50)\n",
    "\n",
    "# Define jobList with initial state |00>\n",
    "for tg in tgList:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job += h1\n",
    "    job += generateIdle(tg, 1)\n",
    "    job += h1\n",
    "    jobListGrd.addJob(job)\n",
    "\n",
    "# Define jobList with initial state |01>    \n",
    "for tg in tgList:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job += x2\n",
    "    job += h1\n",
    "    job += generateIdle(tg, 1)\n",
    "    job += h1\n",
    "    jobListExd.addJob(job)\n",
    "\n",
    "# Run the simulation    \n",
    "stateInit = tensor(basis(level, 0), basis(level, 0), basis(level, 0))\n",
    "resultGrd = runHamiltonian(ham, state0=stateInit, jobList=jobListGrd)\n",
    "resultExd = runHamiltonian(ham, state0=stateInit, jobList=jobListExd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "num0List = []\n",
    "num1List = []\n",
    "\n",
    "# projector |1><1| of qubit 1\n",
    "prj1 = tensor(np.identity(3), basis(level, 1) @ dagger(basis(level, 1)), np.identity(3))\n",
    "\n",
    "# append the result to the list\n",
    "for res0, res1 in zip(resultGrd, resultExd):\n",
    "    psi0, psi1 = res0['state'], res1['state']\n",
    "    rho0, rho1 = psi0 @ dagger(psi0), psi1 @ dagger(psi1)\n",
    "    num0List.append(expect(prj1, rho0))\n",
    "    num1List.append(expect(prj1, rho1))\n",
    "\n",
    "\n",
    "# plot the result\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(tgList, num0List, '.b', label=r'$|00\\rangle$')\n",
    "plt.plot(tgList, num1List, '.r', label=r'$|01\\rangle$')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel('Population of excited state of qubit 1')\n",
    "plt.xlabel\n",
    "plt.title('Ramsey on Q1')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Fit the curve\n",
    "para1Fit, _ = curve_fit(fit, tgList, num0List, [2 * pi * freq / 1.2, 0])\n",
    "para2Fit, _ = curve_fit(fit, tgList, num1List, [2 * pi * freq / 1.2, 0])\n",
    "step = 0.01\n",
    "y1Fit = [fit(x, para1Fit[0], para1Fit[1]) for x in np.arange(tgList[0], tgList[-1], step)]\n",
    "y2Fit = [fit(x, para2Fit[0], para2Fit[1]) for x in np.arange(tgList[0], tgList[-1], step)]\n",
    "\n",
    "# Plot the curve\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(np.arange(tgList[0], tgList[-1], step), y1Fit)\n",
    "plt.plot(np.arange(tgList[0], tgList[-1], step), y2Fit)\n",
    "plt.plot(tgList, num0List, '.b', label=r'$|00\\rangle$')\n",
    "plt.plot(tgList, num1List, '.r', label=r'$|01\\rangle$')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel('Population of excited state of qubit 1')\n",
    "plt.xlabel\n",
    "plt.title('Ramsey on Q1')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "## Calculate the crosstalk strength\n",
    "xiEst = abs(para1Fit[0] - para2Fit[0]) \n",
    "print(f'Coupling strength: {g12 * 1e3 / (2 * pi)} MHz')\n",
    "print(f'ZZ crosstalk strength: {xiEst * 1e3 / (2 * pi)} MHz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从上图可知，因为 $q_1$ 和 $q_2$ 之间耦合强度较小，所以当 $q_2$ 分别处于基态和激发态的时候对应的 $q_1$ 的量子比特频率的差距也不大，说明 $q_1$ 与 $q_2$ 之间的 ZZ 串扰比 $q_0$ 和 $q_1$ 之间的较小。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 总结\n",
    "\n",
    "阅读此关于多比特含噪模拟器的教程后，用户可以通过点击这个链接 [tutorial-multi-qubit-noisy-simulator.ipynb](https://github.com/baidu/Quanlse/blob/main/Tutorial/CN/tutorial-multi-qubit-noisy-simulator-cn.ipynb) 跳转到此 Jupyter Notebook 文档相应的 GitHub 页面，并且运行这个程序。我们鼓励用户使用多比特含噪模拟器探索与本教程不同的其他案例。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 参考文献\n",
    "\\[1\\] [Malekakhlagh, Moein, Easwar Magesan, and David C. McKay. \"First-principles analysis of cross-resonance gate operation.\" *Physical Review A* 102.4 (2020): 042605.](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.102.042605)\n",
    "\n",
    "\\[2\\] [Magesan, Easwar, and Jay M. Gambetta. \"Effective Hamiltonian models of the cross-resonance gate.\" *Physical Review A* 101.5 (2020): 052308.](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.101.052308)\n",
    "\n",
    "\\[3\\] [Ku, Jaseung, et al. \"Suppression of Unwanted ZZ Interactions in a Hybrid Two-Qubit System.\" *Physical review letters* 125.20 (2020): 200504.](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.125.200504)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
